library(specr)
library(lme4)
library(readr)
library(broom.mixed)
library(texreg)
library(dotwhisker)
library(dplyr)

set.seed(42)

source("funs_and_constants.R")

# load in data
epr_yearly_group_clusters_h1 <- read_csv("data/epro_data_h1.csv")
epro_year_group_clusters <- read_csv("data/epro_data_h2.csv")
epro_year_group_clusters_geo <- read_csv("data/epro_geo_data_h3.csv")

# control vector
# groupsize : popshare of group
# incidence flag: Binary  flag indicating ongoing group-level conict
# tek_isrelevant: Binary flag indicating whether at least one of this group’s TEK groups is currently politically relevant
possible_controls <- c("n_orgs", "regaut", "groupsize", "incidence_flag", "status_pwrrank", "any_multiethnic")

############
### H1: resource rents
############
epr_yearly_group_clusters_h1 <- epr_yearly_group_clusters_h1 %>% mutate(across(c(year, countries_gwid, gwgroupid), ~ as.character(.x)))
# convert to character


# random intercept model
lmer_ri <- function(formula, data, ...){
  require(lme4)
  require(broom.mixed)
  formula <- paste(formula, " + (1|countries_gwid) + (1|year) +(1|gwgroupid)")
  lmer(formula, data)
  
}


lmer_ri_no_year <- function(formula, data, ...){
  require(lme4)
  require(broom.mixed)
  formula <- paste(formula, " + (1|countries_gwid) +(1|gwgroupid)")
  lmer(formula, data)
  
}

glmer_ri <- function(formula, data, ...){
  require(lme4)
  require(broom.mixed)
  formula <- paste(formula, "+(1|gwgroupid) + (1|countries_gwid) + (1|year)")
  glmer(formula, data, family = "binomial")
  
}

# check variance of different levels

baseline_model <- lmer(disagreement ~ 1 + (1|gwgroupid) + (1|countries_gwid) + (1|year), epr_yearly_group_clusters_h1)

full_model_h1 <- glmer(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups  + 
                         (1| year) + (1|countries_gwid),
                       data = epr_yearly_group_clusters_h1, family = "binomial")

full_model_h1_max <- lmer(share_disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups +
                            (1| year) + (1 | countries_gwid),
                           data = epr_yearly_group_clusters_h1)

full_model_h1_n_resources <- glmer(disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups  + 
                                     (1| year) + (1|countries_gwid),
                                   data = epr_yearly_group_clusters_h1, family = "binomial")

full_model_h1_n_resources_max <- lmer(share_disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic   + n_tek_groups +
                                        (1| year) + (1|countries_gwid),
                                   data = epr_yearly_group_clusters_h1)


# full_model_h1_no1clusters <- glmer(disagreement ~ resource_and_agriculture + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
#                     +(1|year) + (1|countries_gwid), data = epr_yearly_group_clusters_h1 %>% filter(n_orgs > 1), family = "binomial")

model_list_h1 <- list(full_model_h1, full_model_h1_max, full_model_h1_n_resources, full_model_h1_n_resources_max)

texreg::screenreg(model_list_h1, stars = stars, custom.coef.map = coef_map)

texreg(model_list_h1,
       stars = stars,
       custom.coef.map = coef_map,
       table = FALSE,
       custom.model.names = c("Disagreement (1/0)", "Prop. Disagreement", "Disagreement (1/0)", "Prop. Disagreement"),
       include.proj.stats = FALSE,
       file = "tables/h1_multilevel.tex")

# plot_variance(baseline_model)
# plot_variance(full_model)
# icc_specs(full_model)

#########
# h2
#########

full_model_h2 <- glmer(disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                           +(1|year) + (1|countries_gwid), data = epro_year_group_clusters, family = "binomial")

# n rels as IV
full_model_h2_n_rels <- glmer(disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                       +(1|year) + (1|countries_gwid), data = epro_year_group_clusters, family = "binomial")


full_model_h2_herf <- glmer(disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                       +(1|year) + (1|countries_gwid), data = epro_year_group_clusters, family = "binomial")


# full_model_h2_no1clusters <- glmer(disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
#                        +(1|year) + (1|countries_gwid), data = epro_year_group_clusters %>% filter(n_orgs > 1), family = "binomial")

full_model_h2_max <- lmer(share_disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                       +(1|year) + (1|countries_gwid), data = epro_year_group_clusters)

full_model_h2_max_n_rels <- lmer(share_disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                           +(1|year) + (1|countries_gwid), data = epro_year_group_clusters)

full_model_h2_max_herf <- lmer(share_disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                                  +(1|year) + (1|countries_gwid), data = epro_year_group_clusters)

screenreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf, full_model_h2_max, full_model_h2_max_n_rels, full_model_h2_max_herf),
          stars = stars)

texreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf, full_model_h2_max, full_model_h2_max_n_rels, full_model_h2_max_herf),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h2_multilevel.tex")


#########
# h3
#########

full_model_h3 <- glmer(disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                       + (1|year) + (1|countries_gwid),
                       data = epro_year_group_clusters_geo,
                       family = "binomial")

full_model_h3_max <- lmer(share_disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                       +(1|year) + (1|countries_gwid), data = epro_year_group_clusters_geo)

full_model_h3_bin <- glmer(disagreement ~ multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                           +(1|year) + (1|countries_gwid), data = epro_year_group_clusters_geo, family = "binomial")

full_model_h3_max_bin <- lmer(share_disagreement ~  multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                               +(1|year) + (1|countries_gwid), data = epro_year_group_clusters_geo)

full_model_h3_hhi <- glmer(disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                           +(1|year) + (1|countries_gwid), data = epro_year_group_clusters_geo, family = "binomial")

full_model_h3_max_hhi <- lmer(share_disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                           +(1|year) + (1|countries_gwid), data = epro_year_group_clusters_geo)



full_model_h3_prop_no_outliers <- lmer(agreement_sec_prop ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups
                       +(1|year) + (1|countries_gwid),
                       data = epro_year_group_clusters_geo %>% filter(n_non_intersection_group_polygons < 30))

# full_model_h3_no1clusters <- glmer(disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
#                        +(1|year) + (1|countries_gwid), data = epro_year_group_clusters_geo %>% filter(n_orgs > 1), family = "binomial")

screenreg(list(full_model_h3, full_model_h3_max, full_model_h3_hhi, full_model_h3_max_hhi, full_model_h3_bin, full_model_h3_max_bin), stars = stars)

texreg(list(full_model_h3, full_model_h3_bin, full_model_h3_hhi, full_model_h3_max, full_model_h3_max_bin, full_model_h3_max_hhi),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6),
       table = FALSE,
       file = "tables/h3_multilevel.tex")



######
## commented out as
# specs take long time to run
#######


# make a nice coefficient plot
# coefficient_plot <- plotreg(list(full_model_h1, full_model_h2_no1, full_model_h3_no1),
#         custom.model.names = c("H1", "H2", "H3"),
#         custom.coef.map = c("resource_and_agriculture" = "Diverse Income Streams",
#                             "religious_segments_bin" = "Sectarian Differences",
#                             "n_non_intersection_group_polygons" = "Several Settlement Areas"),
#         custom.note = "",
#         signif.light = "#C5DBE9",
#         signif.medium = "#5A9ECC",
#         signif.dark = "#1C5BA6",
#         insignif.light = "grey",
#         insignif.medium = "grey",
#         insignif.dark = "grey") +
#   theme(text = element_text(size=30))+ ylab("Coefficient Estimate")
# 
# ggsave("plots/coefficient_plot_all.pdf", coefficient_plot, device = "pdf", width = 13.5, height = 8, dpi = 200, units = "in")
# 
# # spec curves for later
# 
# results_h1 <- run_specs(df = epr_yearly_group_clusters_h1,
#                      y = c("disagreement", "share_disagreement"),
#                      x = c("resource_and_agriculture", "none_one_or_both_nats_agri"),
#                      controls = possible_controls,
#                      model = "lmer_ri_no_year", conf.level = 0.9,
#                      all.comb = T
#                      )
# 
# plot_specs(results_h1)
# 
# # higher hhi = more polarization!
# results_h2 <- run_specs(df = epro_year_group_clusters,
#                      y = c("disagreement", "share_disagreement"),
#                      x = c("religious_segments_bin", "hhi_rel", "n_ed_religions"),
#                      controls = possible_controls,
#                      model = "lmer_ri_no_year", conf.level = 0.9,
#                      all.comb = T
# )
# 
# plot_specs(results_h2)
# 
# 
# 
# results_h3 <- run_specs(df = epro_year_group_clusters_geo,
#                         y = c("disagreement", "share_disagreement"),
#                         x = c("n_non_intersection_group_polygons"),
#                         controls = possible_controls,
#                         model = "lmer_ri_no_year", conf.level = 0.9,
#                         all.comb = T
# )
# 
# plot_specs(results_h3)
# 
# # corplot desc
# for_plot <- epro_year_group_clusters_geo %>% select(disagreement, n_non_intersection_group_polygons , n_orgs , regaut , groupsize , incidence_flag , status_pwrrank , any_multiethnic)
# corrplot::corrplot(cor(for_plot, use = "complete.obs"), is.corr = F)






